//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "Unit1.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------
//----------------------------------------------------
//preprocessor directives
#include <stdlib.h>
#include <string>
#include "randomc.h"
#include "userintf.cpp"
#include "mersenne.cpp"
#include "stocc.h"                     // define random library classes
#include "stoc1.cpp"
#include <time.h>
#include <math.h>
#include <fstream>
#include <iostream>
#define xgrid 200   //initial space size
#define ygrid 200
#define xgrid2 6700  // final space size (x dimension expanded)
#define gridsize 10 //note: gridsize must divide perfectly into xgrid, xgrid2, ygrid
#define maxinds2 800000 // max inds in final space
#define PI 3.14159265358979
#define RandomGenerator TRandomMersenne
#define maxdist 20//112.7
#define mutation 0.00001//0.0000329 //mutation rate
#define numsteps 100 //number of steps during dispersal
//-----------------------------------------------------------------------------
//Global structures and variables
struct traits {
  double pdisp;
  double dist;
  double linearity;
  traits();
};
  traits::traits()
  {
    pdisp = 0;
    dist = 0;
    linearity = 0;
  }

struct individual {
  int alive, density, xref, yref;
  double x, y, nextbear, cumdist;
  traits genes;
  individual();
};
  individual::individual(){
    x = 0.0;
    y = 0.0;
    alive = 0;
    density = 0;
    xref = 0;
    yref = 0;
    nextbear = 0.0;
    cumdist = 0.0;
  }
individual inds[maxinds2];
individual temp[maxinds2];
int Xgrid[xgrid2/gridsize*ygrid/gridsize];


const double r=3.0;
//double mutation;
float K;
int popsize;
int32 seed = time(0);
RandomGenerator rg(seed);
StochasticLib1 sto(seed);
/*to avoid runtime error populations shouldn't grow past maxinds*/

//----------------------------------------------------------------------------
//function prototypes
void initialise(void);
void drawgrid(individual jim[], int elements);
void output(individual jim[], int elements);
void output3(int runnum);
void density(individual jim[], int space, int elements);
void reproduction(individual parent[], individual offspring[], int elements); //includes mutation
void dispersal(individual offspring[], int space, int elements);
void cleargrid(individual jo[], int elements);
void transfer(individual parent[], individual offspring[], int elements);
double* means(double array[], int elements);
void jump(void);
//----------------------------------------------------------------------------
//function definitions
void initialise(void){
 int c=0;
 while (c < 0.5*(K*pow(ygrid/gridsize, 2))){   //50% of space K
  inds[c].alive = 1;
  inds[c].x = rg.Random()*xgrid;
  inds[c].y = rg.Random()*ygrid;
  inds[c].genes.pdisp = rg.Random();
  inds[c].genes.dist = rg.Random()*maxdist/numsteps;
  inds[c].genes.linearity = rg.Random(); //proportion of max dist in two moves
  c++;
 }
 popsize=c;
 while (c<maxinds2){
  inds[c].alive = 0;
  inds[c].x = 0;
  inds[c].y = 0;
  inds[c].genes.pdisp = 0;
  inds[c].genes.dist = 0;
  inds[c].genes.linearity = 0;
  c++;
 }
 for (int s=0; s<c; s++){
  Form1->PaintBox1->Canvas->Brush->Color=RGB(0, 0, 0);
  Form1->PaintBox1->Canvas->FillRect(Rect(inds[s].x*2,inds[s].y*2,inds[s].x*2+2,inds[s].y*2+2));
 }
}

void drawgrid(individual jim[], int elements){
 for (int i=0; i<popsize; i++){
  Form1->PaintBox1->Canvas->Brush->Color=RGB(0, 0, 0);
  Form1->PaintBox1->Canvas->FillRect(Rect(jim[i].x*2,jim[i].y*2,jim[i].x*2+2,jim[i].y*2+2));
 }
}

void output(individual jim[], int elements){
 std::ofstream output;
 output.open("testdata.txt", std::ios::app);
 output<<"alive x y xref yref density nextbear cumdist pdisp dist linearity"<<std::endl;
 for (int i=0; i<elements; i++){
  output<<jim[i].alive<<" "<<jim[i].x<<" "<<jim[i].y<<" ";
  output<<jim[i].xref<<" "<<jim[i].yref<<" "<<jim[i].density<<" ";
  output<<jim[i].nextbear<<" "<<jim[i].cumdist<<" "<<jim[i].genes.pdisp;
  output<<" "<<jim[i].genes.dist<<" "<<jim[i].genes.linearity<<std::endl;
 }
 output.close();
}
void density(individual jim[], int space, int elements){
 for (int x=0; x<xgrid2/gridsize*ygrid/gridsize; x++){
  Xgrid[x]=0;
 }
 int yelements = ygrid/gridsize;
 for (int i=0; i<popsize; i++){
  if (jim[i].alive==1){
   jim[i].xref = jim[i].x/gridsize;
   jim[i].yref = jim[i].y/gridsize;
   Xgrid[yelements*jim[i].xref+jim[i].yref]++;
  }
 }
 for (int i=0; i<popsize; i++){
  if (jim[i].alive==1){
   jim[i].density = Xgrid[yelements*jim[i].xref+jim[i].yref];
  }
 }
}
void reproduction(individual parent[], individual offspring[], int elements){
 int tc=0; //counts total number of offspring
 for (int i=0; i<popsize; i++){
 if (parent[i].alive==1){
  double Exp;  //expected num of inds
  int num; //integer number of replacements
  Exp = exp(r*(1-parent[i].density/K));
  Exp = sto.Poisson(Exp);  //introduce stochasticity
  num=Exp;
  if (Exp-Exp/1>=0.5) num++;
  num += tc;
  for (int t=tc; t<num; t++){
   offspring[t].alive=1;
   offspring[t].x = parent[i].x;
   offspring[t].y = parent[i].y;
   offspring[t].nextbear = rg.Random()*360;
   offspring[t].genes.pdisp = parent[i].genes.pdisp;
   offspring[t].genes.dist = parent[i].genes.dist;
   offspring[t].genes.linearity = parent[i].genes.linearity;
   tc++;
   if (sto.Binomial(1, mutation)==1) offspring[t].genes.pdisp = rg.Random();
   if (sto.Binomial(1, mutation)==1) {
    offspring[t].genes.dist = rg.Random()*maxdist/numsteps;
   }
   if (sto.Binomial(1, mutation)==1) offspring[t].genes.linearity = rg.Random();
  }
 }}
 popsize = tc;
 while (tc<maxinds2){       //cleans up offspring array
  offspring[tc].alive = 0;
  tc++;
 }
}
void dispersal(individual offspring[], int space, int elements){
int xsize;
if (space==1) xsize=xgrid;
else xsize=xgrid2;
 for (int t=0; t<popsize; t++){
  double TotDisp;
  if (offspring[t].alive==1){
   TotDisp = numsteps*offspring[t].genes.pdisp*offspring[t].genes.dist*offspring[t].genes.linearity;
   offspring[t].x += TotDisp*cos(offspring[t].nextbear*PI/180);
   offspring[t].y += TotDisp*sin(offspring[t].nextbear*PI/180);
   offspring[t].cumdist = TotDisp/offspring[t].genes.linearity;

   if (offspring[t].x<0) offspring[t].x=offspring[t].x*-1;
   if (offspring[t].x>xsize) offspring[t].x-=2*(offspring[t].x-xsize);
   if (offspring[t].y<0) offspring[t].y=offspring[t].y*-1;
   if (offspring[t].y>ygrid) offspring[t].y-=2*(offspring[t].y-ygrid);

   double PDie;
   if (offspring[t].cumdist==0) PDie = 0;
    else PDie = offspring[t].cumdist/maxdist ;
    //(numsteps/5*1/(1+exp(-0.43711285*log(offspring[t].cumdist)+6.7346219)))+offspring[t].cumdist/maxdist;
   if (PDie>1) PDie=1;
   offspring[t].alive = sto.Binomial(1, 1-PDie);
  }
 }
}
void transfer(individual parent[], individual offspring[], int elements){
 cleargrid(parent, elements);  //function clears drawgrid by ref to inds and counts current population size
 popsize=0;
 for (int i=0; i<elements; i++){
  if (offspring[i].alive==1){
   parent[popsize].alive = offspring[i].alive;
   parent[popsize].x = offspring[i].x;
   parent[popsize].y = offspring[i].y;
   parent[popsize].xref = 0;
   parent[popsize].yref = 0;
   parent[popsize].density = 0;
   parent[popsize].cumdist = 0;
   parent[popsize].nextbear = 0;
   parent[popsize].genes.pdisp = offspring[i].genes.pdisp;
   parent[popsize].genes.dist = offspring[i].genes.dist;
   parent[popsize].genes.linearity = offspring[i].genes.linearity;
   popsize++;
  }
  offspring[i].alive=0;
  offspring[i].x=0;
  offspring[i].y=0;
  offspring[i].xref=0;
  offspring[i].yref=0;
  offspring[i].density=0;
  offspring[i].cumdist=0;
  offspring[i].nextbear=0;
  offspring[i].genes.pdisp=0;
  offspring[i].genes.dist=0;
  offspring[i].genes.linearity=0;
 }
 for (int i=popsize; i<elements; i++){
   parent[i].alive=0;
   parent[i].x=0;
   parent[i].y=0;
   parent[i].xref=0;
   parent[i].yref=0;
   parent[i].density=0;
   parent[i].cumdist=0;
   parent[i].genes.pdisp=0;
   parent[i].genes.dist=0;
   parent[i].genes.linearity=0;
 }
}
void cleargrid(individual jo[], int elements){
 for(int i=0; i<elements; i++){
  Form1->PaintBox1->Canvas->Brush->Color=RGB(255, 255, 255);
  Form1->PaintBox1->Canvas->FillRect(Rect(jo[i].x*2,jo[i].y*2,jo[i].x*2+2,jo[i].y*2+2));
 }
}
double* means(individual jim[], int elements){
 double pdisp, dist, linearity;
 int i=0, c=0;
 static double means[3];
 for (i=0; i<elements; i++){
  if (jim[i].alive==1){
   pdisp += jim[i].genes.pdisp;
   dist += jim[i].genes.dist;
   linearity += jim[i].genes.linearity;
   c++;
  }
 }
 means[0] = pdisp/c;
 means[1] = dist/c;
 means[2] = linearity/c;
 return means;
}
void output3(int runnum){
 int maxX, Xref, origincount, frontcount;
 double dblmaxX, originmean, frontmean;
 dblmaxX=0;
 frontcount=0;
 origincount=0;
 originmean=0;
 frontmean=0;
 for (int c=0; c<popsize; c++){
  if (inds[c].x>dblmaxX) dblmaxX=inds[c].x;
 }
 maxX=dblmaxX/gridsize;
 for (int d=0; d<popsize; d++){
  Xref = inds[d].x/gridsize;
  if (Xref > maxX-5 && Xref <= maxX){
   frontmean += numsteps*inds[d].genes.linearity*inds[d].genes.pdisp*inds[d].genes.dist;
   frontcount++;
  }
  if (Xref < 12 && Xref >= 8){
   originmean += numsteps*inds[d].genes.linearity*inds[d].genes.pdisp*inds[d].genes.dist;
   origincount++;
  }
 }
 frontmean= frontmean/frontcount;
 originmean= originmean/origincount;
 std::ofstream output4;
 output4.open("frontback.txt", std::ios::app);
 output4<<runnum<<" "<<r<<" "<<mutation<<" "<<origincount<<" ";
 output4<<originmean<<" "<<frontcount<<" "<<frontmean<<" "<<maxX<<std::endl;
 output4.close();
}

void jump(void){
 double maxX;
 int maxXref, n;
 maxXref=n=0;
 maxX=0.0;
 for (int i=0; i<popsize; i++){
  if (inds[i].x>maxX) maxX=inds[i].x;
 }
 maxXref = maxX/gridsize;
 n=0;
 int back=5;
 while (n<3000){
  for (int i=0; i<popsize; i++){
   if (inds[i].x/gridsize>(maxXref-back) && inds[i].x/gridsize<=maxXref) n++;
  }
  back++;
 }
 n=0;
 while (n<20){
 for (int i=0; i< popsize; i++){
   if (inds[i].x/gridsize>(maxXref-back) && inds[i].x/gridsize<=maxXref && sto.Binomial(1, 0.01)==1){
    inds[i].x=868.6;
    inds[i].y=190;
    inds[i].density=0.5*K;
    n++;
   }
  }
 }
}


void __fastcall TForm1::Button1Click(TObject *Sender)
{
 /*std::ofstream output4;
 output4.open("frontback.txt", std::ios::app);
 output4<<"run# r mutn origincount originmean frontcount frontmean maxX"<<std::endl;
 output4.close(); */
 std::ofstream output2;
 output2.open("MultiMeansByX.txt", std::ios::app);
 output2<<"initialX finalX GridSize MaxInds2 MaxDist NumSteps K DensDep(local_global) CostInflection CostSlope"<<std::endl;
 output2<<xgrid<<" "<<xgrid2<<" "<<gridsize<<" "<<" "<<maxinds2;
 output2<<" "<<maxdist<<" "<<numsteps<<" "<<K;
 output2<<" local "<<0.75*maxdist<<" "<<0.8<<std::endl<<std::endl;
 output2<<"mutation K gen Xref N pdisp dist linearity run"<<std::endl;
 output2.close();
for (int n=0; n<8; n++){
 switch (n){
  case 0:
   K=5.0;
   break;
  case 1:
   K=10;
   break;
  case 2:
   K=15;
   break;
  case 3:
   K=20;
   break;
  case 4:
   K=25;
   break;
  case 5:
   K=30;
   break;
  case 6:
   K=35;
   break;
  case 7:
   K=40;
   break;
 }
/*for (int m=2; m<3; m++){
 switch (m){
  case 0:
   r=1.05;
   break;
  case 1:
   r=2;
   break;
  case 2:
   r=2.804;
   break;
  case 3:
   r=4;
   break;
 }      */
 for (int z=0; z<10; z++){
  initialise();  //initial values for traits
  for (int c=0; c<100; c++){
   density(inds, 1, maxinds2);   //calculate density in space #1
   /*if (c==99){    //output every ten gens
   output2.open("MultiMeansByX.txt", std::ios::app);
   double pdisp, dist, linearity;
   int N;
    for (int x=0; x<xgrid/gridsize; x++){
    pdisp=0;
    dist=0;
    linearity=0;
    N=0;
     for (int i=0; i<popsize; i++){
      if (inds[i].alive == 1 && inds[i].xref == x){
       pdisp += inds[i].genes.pdisp;
       dist += inds[i].genes.dist;
       linearity += inds[i].genes.linearity;
       N++;
      }
     }
     if (N != 0){
      pdisp = pdisp/N;
      dist = dist/N;
      linearity = linearity/N;
      output2<<mutation<<" "<<r<<" "<<c<<" "<<x<<" "<<N<<" "<<pdisp<<" "<<dist<<" "<<linearity<<" "<<z+1<<std::endl;
     }
    }
    output2.close();
   }  */
   reproduction(inds, temp, maxinds2);  // reproduce (density dependent) + mutation
   dispersal(temp, 1, maxinds2); //disperse with cost linearly related to distance
   drawgrid(temp, maxinds2);
   transfer(inds, temp, maxinds2); // move temp array to final array and clean up
  }

// change population space...
  for (int c=100; c<150/*<=172*/; c++){
   density(inds, 2, maxinds2);
   //if (c==134) jump();
   //if (c%5==0 || c==134 || c==172){
    double pdisp, dist, linearity;
    int N;
     for (int x=0; x<xgrid2/gridsize; x++){
     pdisp=0;
     dist=0;
     linearity=0;
     N=0;
      for (int i=0; i<popsize; i++){
       if (inds[i].alive == 1 && inds[i].xref == x){
        pdisp += inds[i].genes.pdisp;
        dist += inds[i].genes.dist;
        linearity += inds[i].genes.linearity;
        N++;
       }
      }
      if (N!=0){
       pdisp = pdisp/N;
       dist = dist/N;
       linearity = linearity/N;
       output2.open("MultiMeansByX.txt", std::ios::app);
       output2<<mutation<<" "<<K<<" "<<c<<" "<<x<<" "<<N<<" "<<pdisp<<" "<<dist<<" "<<linearity<<" "<<z+1<<std::endl;
       output2.close();
      }
     }
   // }

   reproduction(inds, temp, maxinds2);
   dispersal(temp, 2, maxinds2);
   drawgrid(temp, maxinds2);
   transfer(inds, temp, maxinds2);
  }
//output4();
 }
}//}
}
//---------------------------------------------------------------------------

